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1  Foreword 


This  report  first,  reviews  our  accomplishments  on  the  research  under  AFOSR  Grant  FA9550-08-1- 
0122  in  the  second  year  (March  15,  2009  -  March  14,  2010),  and  then  reviews  the  comprehensive 
accomplishments  over  the  entire  project  period  (Mach  15,  2008  -  March  14,  2010).  Therefore,  this 
report  also  serves  as  the  Final  Report. 

The  high-order  space  time  method  proposed  in  this  research  has  been  termed  as  the  discontinuous 
Galerkin  Cell- Vertex  Scheme  (DG-CVS)  for  its  DG  ingredient,  and  its  alternate  cell  vertex  solution 
updating  strategy. 


2  Accomplishments  in  the  second  year 


In  the  second  year,  the  effort  addressed  the  following  issues: 


1.  Solution  limiting  procedure.  The  limiting  procedure  performs  as  an  alternative  method  to 
limit  (or  remap)  high  order  but.  oscillatory  solutions  across  discontinuities,  thus  maintaining 
the  stability  of  DG-CVS.  In  the  present  approach,  limiting  is  active  only  in  truly  oscillatory 
regions  detected  by  a  reliable  oscillation  indicator.  Since  the  limiter  must  be  conservative, 
the  solution  is  reformulated  in  terms  of  the  cell  average  and  cell-averaged  solution  derivatives, 
and  the  limiter  limits  the  averaged  solution  derivatives  only  while  preserving  the  cell  average. 
The  limiter  ensures  the  solution  satisfies  the  following  two  constraints:  (i)  the  solution  does 
not  exceed  the  maximum  or  minimum  cell  averages  in  the  local  stencil;  and  (ii)  the  solution 
gradient  is  consistent  with  the  local  solution  variation  across  each  edge.  The  limiting  procedure 
is  recast  as  a  quadratic  programming  problem  with  linear  inequality  constraints,  which  can 
be  solved  by  the  active  set  method.  It  is  hoped  that,  through  the  quadratic  programming, 
the  optimum  limiting  factors  can  be  obtained  while  satisfying  the  linear  constraints.  For  high 
order  solutions,  these  constraints  are  formulated  as  the  sufficient  (not  necessary)  conditions. 
The  method  seems  to  work  fine  for  1-D  equations  (systems)  and  extension  to  higher  dimensions 
still  needs  more  work. 
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2.  Employing  the  quadruple  precision  arithmetic  via  MPACK  |l).  The  purpose  of  doing  this  is  to 
verify  that  the  DG-CVS  is  truly  arbitrarily  high-order  accurate  as  long  as  the  linear  system 
resulted  from  discretization  can  be  accurately  solved.  In  the  current  implementation,  Taylor 
polynomials  are  used  as  the  basis  functions  and  Taylor  polynomials  are  notorious  in  generating 
severely  ill-conditioned  systems.  Therefore,  double-precision  arithmetic  cannot  reliably  solve 
the  system  when  the  degree  (p)  of  basis  polynomials  is  high  (p  >  4).  The  quadruple  precision 
arithmetic  provided  by  MPACK  is  able  to  solve  highly  ill-conditioned  systems  and  verifies 
that  optimal  convergence  rate,  p  +  1,  can  be  readied  for  p  >  4. 

3.  Improving  the  efficiency  of  the  present  high-order  space-time  method  by  implementing  it  on 
the  overset  Cartesian/ quadrilateral  grids.  Though  the  DG-CVS  is  designed  for  arbitrarily 
unstructured  grids,  rectangular /quadrilateral  grids  are  superior  to  triangular  grids  in  terms 
of  efficiency.  The  body-fitted  quadrilateral  mesh  is  used  to  wrap  the  object  and  the  Cartesian 
grid  serves  as  the  background  mesh.  The  contact  boundary  extracted  from  the  quadrilateral 
mesh  is  used  to  determine  the  overlapping  region  between  the  two  meshes.  The  overlapping 
region  is  kept  minimum  but  large  enough  to  transfer  solutions  between  the  two  meshes.  The 
Cartesian  mesh  around  the  overlapping  region  is  geometrically  refined  to  match  the  resolution 
of  the  corresponding  quadrilateral  mesh.  The  portion  of  the  Cartesian  mesh  covered  by  the 
quadrilateral  mesh  and  the  object  and  beyond  the  the  overlapping  region  is  discarded.  The 
inter-mesh  solution  transfer  strategy  is  crucial  in  ensuring  the  global  flux  conservation.  The 
donor  cell  is  the  cell  on  the  counterpart  mesh  where  the  receptor  cell  is  located.  A  simple  test 
of  a  subsonic  flow  around  a  circular  cylinder  demonstrates  that  the  current  implementation  is 
quite  successful. 

4.  Dissemination  via  conference  and  publications.  The  limiting  procedure  described  above  has 
been  presented  and  published  (AIAA  Paper  2009-3983)  in  the  19th  AIAA  Computational 
Fluid  Dynamics  Conference,  June,  2009,  San  Antonio,  TX.  The  results  of  the  implementation 
on  overset  Cartesiau/quadrilateral  grids  have  been  presented  and  published  (AIAA  Paper 
2010-0544)  in  the  48th  AIAA  Aerospace  Science  Meeting,  January  2010  in  Orlando,  Florida. 
A  revised  comprehensive  paper  is  to  appear  in  Communications  in  Computational  Physics. 


3  Accomplishments  in  the  entire  grant  period 


Over  the  entire  grant  period,  the  present  high-order  space  time  method,  which  is  termed  as  the 
discontinuous  Galerkin  cell- vertex  scheme  (DG-CVS),  has  been  verified  and  developed  for  both 
scalar  advection  equations  and  compressible  Euler  equations.  The  work  so  far  shows  that  DG-CVS 
is  a  promising  method  for  hyperbolic  conservation  laws.  The  most  distinct  features  of  the  DG-CVS 
can  be  sununarized  as: 


•  locally  and  globally  space-time  flux  conservative. 

•  alternate  solution  updating  at  the  cell  level  and  the  vertex  level  within  each  physical  time 
step. 

•  Riemann-solver-free  for  advection  problem. 

•  higli-order  accurate  in  both  space  and  time. 

•  liiglily  compact  regardless  of  the  order. 
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•  suitable  for  arbitrarily  unstructured  meshes. 

•  simple  boundary  condition  treatment. 


A  comprehensive  list  of  accomplishments  is  given  below: 


1.  Alternate  cell  vertex  solution  updating  strategy.  Inspired  by  the  space-time  conservation  el¬ 
ement/solution  element  (CE/SE)  method  |2),  space-time  flux  conservation  is  enforced  over 
staggered  space-time  conservation  elements.  In  this  research  universal  definitions  of  CEs  and 
SIDs  using  the  cell-vertex  structure  of  the  underlying  spatial  mesh  are  utilized  such  that  the 
method  is  suitable  for  arbitrarily  unstructured  meshes.  The  solution  within  each  physical 
time  step  is  updated  alternately  at  the  cell  level  and  the  vertex  level.  Tliis  strategy  cir¬ 
cumvents  the  need  of  using  Riemann  solvers  to  provide  the  inter-cell  fluxes  and  leads  to  a 
Riemann -solver-free  solver. 


2.  Arbitrarily  high-order  accuracy  in  both  space  and  time.  The  high  order  of  accuracy  is  achieved 
by  employing  high-order  polynomials  basis  functions  inside  each  SE,  as  in  the  discontinuous 
Galerkin  (DG)  method  |3|.  What  is  different  from  the  classic  DG  method  is  that  no  any  type 
of  Riemann  solvers  is  used  for  the  inter-cell  flux  thanks  to  the  cell  vertex  solution  updating 
strategy  explained  above.  This  DG  ingredient  also  makes  the  current  method  deviate  from 
the  CE/SE  method  by  providing  high  order  accuracy.  For  the  DG  ingredient  and  the  cell- 
vertex  solution  updating  strategy,  the  new  scheme  here  is  termed  as  the  discontinuous  Galerkin 
cell- vertex  scheme  (DG-CVS). 


3.  Formulation.  Considering  the  following  one  dimension  linear  scalar  advection  equation 


du{xA)  ,  df(u) 


dt 


dx 


=  0 


(1) 


where  u  is  the  advected  quantity  and  /  is  the  flux.  The  approximate  solution  uh  is  sought 
within  each  space-time  solution  element  (SE),  denoted  as  K.  When  restricted  to  the  SE,  uh 
belongs  to  the  finite  dimensional  space  U(K)  such  that 


«*(*.«)  =  (2) 

3= 1 


where  are  some  type  of  polynomial  basis  functions,  are  the  unknowns  to  be 

determined  and  N  is  the  number  of  basis  functions  depending  on  the  degree  of  the  polynomial 
function. 

According  to  Galerkin  orthogonality,  one  can  multiply  the  governing  equation  with  each  of 
the  basis  functions  {0,};^  to  obtain 

ir+i ^■)dn  =  "  1 °“=v-,N  (3) 

where  Q  is  the  conservation  element  (CE)  corresponding  to  the  solution  element  K. 
Integrating  the  resulting  weak  form  by  parts  yields 


(4) 
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Figure  1:  Illustration  of  space-time  flux  conservation  in  a  conservation  element  at  the  cell 
level. 


where 

nx,nt)  (5) 

is  the  space-time  flux  normal  to  the  boundary  of  the  space-time  CE.  n  =  (nx,  n,)  is  the  outward 
unit  normal  of  the  CE  boundary.  T  =  dfl  is  the  boundary  of  the  CE  under  consideration. 


The  cell  level  CE  shown  in  Fig.  1  is  taken  as  a  specific  example  here.  As  shown  in  Fig.  1, 
divide  T  into  five  sections  rt,  T2,  r3,  T.,  and  T5  where  T j  belongs  to  the  SE  associated  with 
(m  +  2.  n  +  j)  where  the  solution  is  being  sought,  T2  and  r3  the  SE  associated  with  (m,n) 
and  T,  and  Ts  the  SE  associated  with  (m  +  1,  n).  Note  that  the  solutions  at  nodes  (m,  n)  and 
(m  +  1,  n)  are  known  since  they  are  at  the  previous  time  level.  Considering  the  outward  unit 
normal  of  each  boundary,  Eq.  (4)  becomes 


which  leads  to  a  linear  or  nonlinear  equation  system  (depending  on  whether  /  is  a  linear  or 
nonlinear  function  of  u)  that  can  be  solved  for  the  unknown  u^+1^2  n+1y2  at  the  space-time 
node  (m  +  1/2,  n  +  1/2). 


4.  Choice  of  basis  polynomials.  In  this  project  ,  Taylor  polynomials  are  chosen  as  the  basis  func¬ 
tions  for  tliree  reasons.  First,  the  Taylor  polynomial  has  no  restrictions  on  the  geometric 
shape  of  the  conservation  element  (SE).  The  SE  can  be  of  arbitrary  shape  which  is  typically 
polygonal  cylinder  on  spatially  unstructured  meshes  .  By  contrast,  Lagrange  polynomials  are 
only  well  defined  on  simplicial  elements  or  elements  allowing  tensor  products,  and  Chebysliev 
polynomials  and  Legendre  polynomials  are  only  well  defined  on  elements  allowing  tensor  prod¬ 
ucts. 

Second,  with  Taylor  polynomials,  the  basis  functions  are  polynomials  with  respect  to  the 
physical  coordinates  Ax  and  At.  Other  polynomials  usually  define  the  basis  function  using 
reference  coordinates  (e.g.,  £  and  //).  The  numerical  integration  involves  the  derivatives  of 
the  basis  function  with  respect  to  the  physical  coordinates.  With  Taylor  polynomials,  the 
derivatives  can  be  directly  obtained  by  taking  the  derivative  of  the  basis  functions  with  re¬ 
spect  to  Ax  and  At  without  resorting  to  the  chain  rule  as  is  required  when  other  polynomials 
are  employed. 

Third,  the  derivatives  of  Taylor  polynomials  are  a  lower  order  subset  of  the  original  Taylor 
polynomials.  Tliis  allows  efficient  implementation  when  evaluating  integrals  involving  prod¬ 
ucts  of  polynomials. 
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The  above  discussions  of  the  advantages  of  Taylor  polynomials  over  other  polynomials  are 
completely  based  on  the  implementation  point  of  view.  One  should  keep  in  mind  that  high 
degree  Taylor  polynomials  are  notorious  in  generating  severely  ill-conditioned  systems.  It  is 
well  known  that  the  number  of  reliable  digits  of  the  solution  of  linear  systems  is  related  to 
the  conditioning  of  the  system.  As  a  result,  high  degree  Taylor  polynomials  are  expected  to 
yield  inaccurate  solutions  due  to  the  ill  conditioning  of  the  system.  Therefore,  choosing  Taylor 
polynomials  is  not  wise  from  the  accuracy  point  of  view.  However,  the  goal  of  this  research 
is  to  demonstrate  the  efficacy  of  the  DG-CVS  idea,  so  Taylor  polynomials  are  employed  for 
quick  implementation. 

5.  Insensitivity  to  Courant  number  when  using  high-order  basis  polynomials.  It  is  well  known  that 
some  space-time  and  fully  discrete  methods  exhibits  excess  dissipation  when  the  time  step  is 
vanishingly  small.  The  present  study  shows  that  the  present  DG-CVS  pi  case  exhibits  similar 
phenomenon,  namely,  when  the  time  step  is  vanishingly  small,  the  accuracy  is  degraded. 
However,  for  higher  degrees,  the  DG-CVS  is  less  and  less  sensitive  to  small  Courant  numbers. 
Actually,  in  the  case  of  p4,  a  slightly  liigher  accuracy  is  obtained  when  a  smaller  time  step  is 
used. 

6.  Quadrature-free  implementation.  The  DG-CVS  formulation  involves  both  surface  and  volume 
integrals.  If  the  underlying  spatial  mesh  is  unstructured,  the  vertex-level  CEs  are  general 
polygonal  cylinders  containing  general  polygonal  bases  and  quadrilateral  side  faces  where  the 
Gaussian  quadrature  rule  cannot  be  directly  applied.  In  addition,  even  for  purely  rectan¬ 
gular  meshes  where  the  Gaussian  quadrature  rule  can  be  applied  to  all  integrals,  integra¬ 
tion  by  quadrature  rule  is  prohibitively  expensive,  especially  in  high  dimensions.  Therefore, 
quadrat ure-free  implementation  is  sought  in  the  research  to  improve  the  efficiency  of  the 
method.  We  can  combine  the  divergence  theorem  and  indefinite  integration  to  transform  a 
volume  integral  to  surface  integrals  and  further  to  line  integrals.  Finally,  analytical  formulae 
are  available  to  evaluate  the  line  integrals.  Numerical  experiments  show  that  the  quadrature- 
free  implementation  drastically  increases  the  efficiency  by  a  factor  of  60  for  the  case  of  p  =  4. 

7.  Boundary  condition  treatment.  The  space-time  formulation  makes  it  simple  to  implement  the 
Dirichlet  and  Neumann  boundary  conditions.  The  boundary  conditions  are  accounted  for 
by  modifying  the  resulting  linear  equation  system  from  discretization.  Both  the  left  hand 
side  matrix  and  right  hand  side  vector  are  modified  consistently.  For  the  outflow  boundary 
condition,  nothing  needs  to  be  modified  as  long  as  the  matrix  is  formed  by  taking  into  account 
the  contribution  from  the  outer  side  face  of  the  boundary  CE. 

8.  Accuracy  for  linear  scalar  advection  equations  (one  dimension  or  two  dimensions).  In  linear 
advection  equations,  the  flux  satisfies  the  same  polynomial  distribution  as  the  solution  itself 
scaled  by  a  constant  advection  speed.  If  the  solution  inside  the  conservation  element  is  assumed 
to  satisfy  a  polynomial  of  degree  p,  then  the  order  of  accuracy  should  bep-t-  1.  This  has  been 
confirmed  in  the  grid  convergence  study  for  both  1-D  and  2-D  problems.  The  convergence 
rates  are  truly  optimal  p  +  1  as  long  as  the  equation  system  resulted  from  discretization  can 
be  solved  accurately.  When  p  >  4,  the  system  becomes  highly  ill-conditioned  due  to  the  use 
of  Taylor  polynomial  bases.  Therefore,  quadruple  precision  arithmetic  is  required  to  obtain 
the  optimal  convergence  rate  when  p  is  high. 

9.  Accuracy  for  nonlinear  scalar  advection  equations  (one  dimension  or  two  dimensions).  In 
nonlinear  advection  equations  (e.g.,  inviscid  Burgers  equation),  the  flux  is  a  nonlinear  function 
of  the  solution.  Therefore,  the  flux  satisfies  a  different  polynomial  distribution  than  the 
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1 1 error  |  oo  under  different  Couraut  number 


Conrant  number  a  =  \a\Sl/6x 

Pi 

p2 

p3 

p4 

0.2 

3.31E-02 

1.42E-04 

2.62E-05 

4.12E-07 

0.02 

1.82E-01 

3.18E-04 

2.63E-05 

2.42E-07 

0.002 

8.20E-01 

3.14E-03 

2.83E-05 

2.08E-07 

Table  1:  Sensitivity  to  the  Couraut  number  of  the  present  method. 


solution  itself.  Thanks  to  the  quadrature- free  implementation,  the  integral  involving  the  flux 
can  be  evaluated  exactly.  The  grid  convergence  study  shows  that  the  optimal  convergence 
rate  of  the  solution  can  also  be  reached  for  inviscid  1-D  and  2-D  Burgers  equations. 

10.  Accuracy  on  unstructured  triangular  meshes.  The  grid  convergence  study  has  also  been  con¬ 
ducted  on  unstructured  triangular  meshes.  Again,  optimal  convergence  rate  can  be  reached 
for  each  p. 

11.  Extension  to  compressible  Euler  equations  (one  dimension  or  two  dimensions).  The  DG-CVS 
idea  can  be  directly  applied  to  more  complex  hyperbolic  compressible  Euler  equations.  In  our 
implementation,  following  Lowrie  et  al.  (4),  we  choose  the  working  variables  to  be 

Q  =  \Sp,  Vp",  Vi*-',  VpH)t  (7) 

where  p,  u,  v,  II  are  density,  x-  and  y- velocity  components,  and  specific  total  enthalpy,  re¬ 
spectively.  By  choosing  such  working  variables,  all  components  of  the  conservative  state  vector 
and  flux  vectors  can  be  expressed  as  the  double  product  between  working  variables,  which  is 
analogous  to  the  inviscid  Burger’s  equation. 

12.  Solution  limiting  procedure.  This  has  been  explained  in  the  section  of  “Accomplishments  in 
the  second  year”. 

13.  Quadruple  precision  implementation.  This  has  been  explained  in  the  section  of  “Accomplish¬ 
ments  in  the  second  year". 

14.  Implementation  on  overset  Cartesian/ quadrilateral  meshes.  This  has  been  explained  in  the 
section  of  “Accomplishments  in  the  second  year”. 


4  Numerical  illustrations 


This  section  presents  some  numerical  results  to  illustrate  the  findings  fisted  in  sections  2  and  3. 

Table  1  lists  the  oc-errors  of  pi,  p2,  p3  and  p4  cases  when  various  Couraut  numbers  are  used  for 
a  linear  advection  equation. 

Figure  2  shows  the  convergence  order  of  the  DG-CVS.  It  can  be  seen  that  the  optimal  rate,  p+ 1, 
can  be  reached  for  all  p  (degree  of  the  basis  polynomial).  Note  that  the  p5  and  p6  cases  are  run 
using  the  quadruple  precision  arithmetic  due  to  the  ill- conditioning  caused  by  Taylor  polynomials. 
Figure  shows  similar  phenomenon  for  nonlinear  inviscid  Burgers  equation.  Figures  4  and  5  show 
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Table  2:  Order  of  accuracy  for  the  1-D  linear  sinusoidal  wave  advection  at  t 
precision. 

p  Tic  h  error  order  loo  error  order 

_ 4  10  2.80E-07  I  4.12E-07  - 

20  8.18E-09  5.096  1.25E-08  5.037 

40  2.46E-10  5.057  3.81E-10  5.041 

80  7.48E-12  5.039  1.17E-11  5.029 

160  2.32E-13  5.013  3.63E-13  5.008 

5  I  10  1.31E-08  "  2.19E-08  I 

20  2.14E-10  5.933  3.53E-10  5.960 

40  3.42E-12  5.967  5.51E-12  6.001 

80  5.39E-14  5.986  8.58E-14  6.004 

160  8.44E-16  5.999  1.33E-15  6.007 

~~ 6  10-  2.00E-11  -  2.95E-11  " 

20  1.43E-13  7.129  2.16E-13  7.094 

40  1.06E-15  7.072  1.64E-15  7.039 

80  8.18E-18  7.022  1.28E-17  7.009 

160  6.38E-20  7.003  9.98E-20  6.997 


=  1.0.  Quadruple 


Table  3:  Order  of  accuracy  for  the  1-D  iuviscid  Burgers  equation  at  t  =  0.12.  Quadruple  precision. 


It  error 


5.57E-06 

1.79E-07 

8.09E-09 

2.47E-10 

7.60E-12 


1.44E-06 
1 .84E-08  6.296 

5.33E-10  5.108 

9.86E-12  5.756 

1.67E-13  5.882 


4.98E-07 
6.82E-09  6.188 

5.87E-11  6.862 

5.15E-13  6.832 

4.61E-15  6.803 


6.93E-06 
1.16E-07 
6.81E-09 
1.18E-10  5.848 

1.89E-12  5.968 


1.80E-06 
3.23E-08  5.799 

3.55E-10  6.510 

4.46E-12  6.312 

4.63E-14  6.592 


160 


Table  4:  Order  of  accuracy  for  the  2-D  linear  sinusoidal  wave  advection  at  t  =  1.0  on  unstructured 
triangular  meshes. 


p 

mesh  size  (/*) 

li  error 

order 

loo  error 

order 

1 

0.2 

5.21E-02 

- 

1.07E-01 

- 

0.1 

9.14E-03 

2.51 

2.01E-02 

2.41 

0.05 

1.70E-03 

2.43 

3.97E-03 

2.34 

0.025 

3.48E-04 

2.29 

8.37E-04 

2.25 

0.0125 

7.71E-05 

2.17 

1.91E-04 

2.13 

2 

0.2 

4.68E-04 

- 

1.90E-03 

- 

0.1 

4.06E-05 

3.53 

2.94E-04 

2.69 

0.05 

3.88E-06 

3.39 

3.60E-05 

3.03 

0.025 

4.09E-07 

3.25 

4.46E-06 

3.01 

0.0125 

4.60E-08 

3.15 

5.60E-07 

2.99 

3 

0.2 

7.05E-05 

- 

2.77E-04 

- 

0.1 

4.48E-06 

3.98 

1.76E-05 

3.98 

0.05 

2.84E-07 

3.98 

1.11E-06 

3.99 

0.025 

1.80E-08 

3.98 

6.96E-08 

4.00 

4 

0.2 

1.91E-06 

- 

8.61E-06 

- 

0.1 

5.79E-08 

5.04 

2.80E-07 

4.94 

0.05 

1.75E-09 

5.05 

9.41E-09 

4.90 

0.025 

5.39E-11 

5.02 

3.18E-10 

4.89 

Table  5:  Order  of  accuracy  for  the  2-D  inviscid  Burgers  equation  at  t  =  0.075  on  unstructured 
triangular  meshes. 


V 

■□rah  size  (h) 

"iter 

l\  error 

order 

loo  error 

order 

1 

0.2 

3.39 

2.12E-02 

- 

9.92E-02 

- 

0.1 

3.53 

5.58E-03 

1.93 

2.79E-02 

1.83 

0.05 

3.18 

1.40E-03 

1.99 

7.26E-03 

1.94 

0.025 

2.95 

3.46E-04 

2.02 

1.75E-03 

2.05 

0.0125 

2.96 

8.59E-05 

2.01 

4.26E-04 

2.04 

2 

0.2 

3.38 

7.50E-04 

- 

1.14E-02 

- 

0.1 

3.55 

1.03E-04 

2.86 

1.73E-03 

2.72 

0.05 

3.17 

1.19E-05 

3.11 

2.60E-04 

2.73 

0.025 

2.94 

1.21E-06 

3.30 

4.21E-05 

2.63 

0.0125 

2.97 

1.28E-07 

3.24 

7.14E-06 

2.56 

3 

0.2 

3.52 

4.35E-04 

- 

6.34E-03 

- 

0.1 

3.64 

3.47E-05 

3.65 

3.85  E-04 

4.04 

0.05 

3.29 

2.42E-06 

3.84 

4.00E-05 

3.27 

0.025 

2.95 

1.52E-07 

3.99 

2.73E-06 

3.87 

0.0125 

2.97 

9.55E-09 

3.99 

1.67E-07 

4.03 

4 

0.2 

3.58 

6.10E-05 

- 

7.15E-04 

- 

0.1 

3.62 

2.71E-06 

4.49 

3.82E-05 

4.23 

0.05 

3.26 

8.18E-08 

5.05 

1.77E-06 

4.43 

0.025 

2.94 

2.36E-09 

5.12 

8.07E-08 

4.46 
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to  (b) 

Figure  2:  Illustration  of  solution  limiting. 


Table  6:  Comparison  of  the  time  cost  in  integration  using  the  Gaussian  quadrature  rule  and 
quadrature-free  approach. 


Scaled  time  cast. 


1  P1 

p2 

| 

| 

Compute  j \ItJ  contributions 
from  the  top  face  of  CEs 

Gaussian  quadrature 

|  0.90 

1.26 

2.08  | 

|  4.34 

Quadrature-free 

j  1 

1 

1  | 

j  1 

Compute  iVf(J  contributions 

from  the  volume  of  CEs 

Gaussian  quadrature 

|  1.78 

5.64 

24.48  1 

1  60.38 

Quadrature-free 

j  1 

1 

1  j 

j  1 

the  convergence  orders  of  2-D  linear  and  nonlinear  advectiou  equations  on  unstructured  triangular 
meshes. 

Figure  2  illustrates  the  performance  of  solution  limiting. 

Figure  3  shows  the  boundary  condition  treatment  of  inflow  and  outflow  boundary  conditions. 

Table  6  illustrates  the  cost  saving  using  the  quadrature-free  implementation  compared  to  the 
Gaussian  quadrature  integration. 

Figures  4,  5  and  6  present  some  results  for  compressible  Euler  equations  including  the  1-D  shock 
tube  problem,  2-D  iseutropic  vortex  advectiou  problem  and  the  forward- facing  step  problem. 

Figure  7  show's  the  simulated  pressure  field  on  an  overset  Cartesiau/quadrilateral  mesh  of  a 
subsonic  flow  around  a  circular  cylinder. 


9 


(a)  t  =  0.0 

(b)  t  =  0.2 

(c)  t  =  0.4 

(d)  t  =  0.6 

(e)  t  =  0.8 

(f)  t  =  1.0 

Figure  3:  Advection  of  a  doubly  raised  cosine  surface.  p2  solution  on  40  x  40  rectangular  mesh. 
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0.9 
0.8 
0.7 
a.  0.6 
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0.3 
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0.1 
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Figure  4:  Density  distribution  of  the  1-D  shock  tube  problem  at  t  =  0.2  when  various  degrees  of 
basis  polynomials  employed.  Without  limiter. 
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Figure  5:  Advection  of  isentropic  vortex  ou  a  20  x  20  rectangular  mesh.  Comparison  of  density 
accuracy  between  pi  —  p4  cases.  The  close-up  view  is  ou  the  right. 


Figure  6:  Solution  of  a  supersonic  (M  =  3)  flow  through  a  channel  with  forward- facing  step  at 
I.  =  4.0.  Without  limiter.  Left:  pi  density  field.  Right:  residual  of  the  continuity  equation. 


Figure  7:  Steady  pressure  field  the  inviscid  flow  (M  =  0.1)  passing  around  a  circular  cylinder. 
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